function [ vh ] = vcycle( vh, fh, current)
%VCYCLE Summary of this function goes here
%   Detailed explanation goes here

global A;
global P;
global v1;
global v2;
global k;
global omega;
%global dimm; % for debug

% Step 1: Relax v1 times on Au=f
for i = 1 : v1
    vh=jacobi(vh, fh, current);
end

% Step 2: If current iteration is not the last one:
if current < k
    f2h = P{current}'*(fh-A{current}*vh);
    v2h = zeros(length(f2h), 1);
    v2h = vcycle(v2h, f2h, current+1);
    % Step 3: Correct vh:
    vh = vh + P{current}*v2h; 
else
    omega_saved = omega;
    omega = 1;
    % Do many relaxations on coarsest grid to solve well
    for i = 1:10
        vh=jacobi(vh, fh, current);
    end
    omega = omega_saved;
end

% Step 4: Relax v2 times on Au=f
for i = 1 : v2
    vh=jacobi(vh, fh, current);
end

end
